#Set working directory to project directory and load packages

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())

library(knitr)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.1     ✔ tibble    3.1.8
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## 
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.18.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## 
## The following object is masked from 'package:survival':
## 
##     kidney
## 
## The following object is masked from 'package:stats':
## 
##     ar
library(parallel)
library(ggridges)
library(loo)
## This is loo version 2.5.1
## - Online documentation and vignettes at mc-stan.org/loo
## - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
library(bayesplot)
## This is bayesplot version 1.10.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting

Load data and set instructions for sampler

ncrmp <- read_csv("data/outputs/ncrmp_simple.csv")
## New names:
## Rows: 1637 Columns: 49
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (5): SampleID, REGION, SUB_REGION_NAME, rugosityMethod, SEASON dbl (44): ...1,
## totalBiomass, coralAssociated, herb, otherFish, YEAR, MONTH,...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
logit <- function(mu){
  logit.mu <- log(mu/(1-mu))
}

#-Sampler instructions-#
n_cores <- detectCores() 
n_chains <- 2 
n_iter <- 4000
n_warmup <- 2000

Build individual models

SCTLD-susceptible coral

################################################################################
#-Build individual models-######################################################
################################################################################

#-Susceptible coral-############################################################

suscCoral_mod <- bf(propSusc ~  1 + 
                      (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                    phi ~ 1 + 
                      (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                    zi ~ 1 + 
                      (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                    family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

get_prior(suscCoral_mod, data = ncrmp)
##                 prior     class         coef           group resp dpar nlpar lb
##                (flat)         b                                                
##                (flat)         b sMAX_DEPTH_1                                   
##  student_t(3, 0, 2.5) Intercept                                                
##  student_t(3, 0, 2.5)        sd                                               0
##  student_t(3, 0, 2.5)        sd              SUB_REGION_NAME                  0
##  student_t(3, 0, 2.5)        sd    Intercept SUB_REGION_NAME                  0
##  student_t(3, 0, 2.5)        sd                         YEAR                  0
##  student_t(3, 0, 2.5)        sd    Intercept            YEAR                  0
##  student_t(3, 0, 2.5)       sds                                               0
##  student_t(3, 0, 2.5)       sds s(MAX_DEPTH)                                  0
##                (flat)         b                                    phi         
##                (flat)         b sMAX_DEPTH_1                       phi         
##  student_t(3, 0, 2.5) Intercept                                    phi         
##  student_t(3, 0, 2.5)        sd                                    phi        0
##  student_t(3, 0, 2.5)        sd              SUB_REGION_NAME       phi        0
##  student_t(3, 0, 2.5)        sd    Intercept SUB_REGION_NAME       phi        0
##  student_t(3, 0, 2.5)        sd                         YEAR       phi        0
##  student_t(3, 0, 2.5)        sd    Intercept            YEAR       phi        0
##  student_t(3, 0, 2.5)       sds                                    phi        0
##  student_t(3, 0, 2.5)       sds s(MAX_DEPTH)                       phi        0
##                (flat)         b                                     zi         
##                (flat)         b sMAX_DEPTH_1                        zi         
##        logistic(0, 1) Intercept                                     zi         
##  student_t(3, 0, 2.5)        sd                                     zi        0
##  student_t(3, 0, 2.5)        sd              SUB_REGION_NAME        zi        0
##  student_t(3, 0, 2.5)        sd    Intercept SUB_REGION_NAME        zi        0
##  student_t(3, 0, 2.5)        sd                         YEAR        zi        0
##  student_t(3, 0, 2.5)        sd    Intercept            YEAR        zi        0
##  student_t(3, 0, 2.5)       sds                                     zi        0
##  student_t(3, 0, 2.5)       sds s(MAX_DEPTH)                        zi        0
##  ub       source
##          default
##     (vectorized)
##          default
##          default
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##          default
##     (vectorized)
##          default
##     (vectorized)
##          default
##          default
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##          default
##     (vectorized)
##          default
##     (vectorized)
##          default
##          default
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##          default
##     (vectorized)
propSusc_mean <- ncrmp %>%
  dplyr::group_by(SUB_REGION_NAME) %>%
  dplyr::summarize(mean = mean(propSusc)) %>%
  mutate(logitmean = logit(mean))

print(paste("The logit-transformed propSusc intercept is", round(propSusc_mean[1,3], 2)))
## [1] "The logit-transformed propSusc intercept is -2.06"
# "The logit-transformed propSusc intercept is -2.06"

susc_prior <- c(prior(normal(-2.1, 1), class = "Intercept"),
                prior(normal(0, 1), class = "b"),
                prior(normal(0, 2), class = "sds"),
                prior(exponential(2), class = "sd"),
                prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
                prior(exponential(2), class = "sd", dpar = "zi"),
                prior(normal(0, 2), class = "sds", dpar = "zi"))


susc_prior$resp <- "propSusc"

SCTLD-resistant coral

#-Resistant coral-##############################################################

resistantCoral_mod <- bf(propResistant ~  1 + 
                           (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                         phi ~ 1 + 
                           (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                         zi ~ 1 + 
                           (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                         family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

get_prior(resistantCoral_mod, data = ncrmp) 
##                 prior     class         coef           group resp dpar nlpar lb
##                (flat)         b                                                
##                (flat)         b sMAX_DEPTH_1                                   
##  student_t(3, 0, 2.5) Intercept                                                
##  student_t(3, 0, 2.5)        sd                                               0
##  student_t(3, 0, 2.5)        sd              SUB_REGION_NAME                  0
##  student_t(3, 0, 2.5)        sd    Intercept SUB_REGION_NAME                  0
##  student_t(3, 0, 2.5)        sd                         YEAR                  0
##  student_t(3, 0, 2.5)        sd    Intercept            YEAR                  0
##  student_t(3, 0, 2.5)       sds                                               0
##  student_t(3, 0, 2.5)       sds s(MAX_DEPTH)                                  0
##                (flat)         b                                    phi         
##                (flat)         b sMAX_DEPTH_1                       phi         
##  student_t(3, 0, 2.5) Intercept                                    phi         
##  student_t(3, 0, 2.5)        sd                                    phi        0
##  student_t(3, 0, 2.5)        sd              SUB_REGION_NAME       phi        0
##  student_t(3, 0, 2.5)        sd    Intercept SUB_REGION_NAME       phi        0
##  student_t(3, 0, 2.5)        sd                         YEAR       phi        0
##  student_t(3, 0, 2.5)        sd    Intercept            YEAR       phi        0
##  student_t(3, 0, 2.5)       sds                                    phi        0
##  student_t(3, 0, 2.5)       sds s(MAX_DEPTH)                       phi        0
##                (flat)         b                                     zi         
##                (flat)         b sMAX_DEPTH_1                        zi         
##        logistic(0, 1) Intercept                                     zi         
##  student_t(3, 0, 2.5)        sd                                     zi        0
##  student_t(3, 0, 2.5)        sd              SUB_REGION_NAME        zi        0
##  student_t(3, 0, 2.5)        sd    Intercept SUB_REGION_NAME        zi        0
##  student_t(3, 0, 2.5)        sd                         YEAR        zi        0
##  student_t(3, 0, 2.5)        sd    Intercept            YEAR        zi        0
##  student_t(3, 0, 2.5)       sds                                     zi        0
##  student_t(3, 0, 2.5)       sds s(MAX_DEPTH)                        zi        0
##  ub       source
##          default
##     (vectorized)
##          default
##          default
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##          default
##     (vectorized)
##          default
##     (vectorized)
##          default
##          default
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##          default
##     (vectorized)
##          default
##     (vectorized)
##          default
##          default
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##          default
##     (vectorized)
propResistant_mean <- ncrmp %>%
  dplyr::group_by(SUB_REGION_NAME) %>%
  dplyr::summarize(mean = mean(propResistant)) %>%
  mutate(logitmean = logit(mean))

print(paste("The logit-transformed propResistant intercept is", round(propResistant_mean[1,3], 2)))
## [1] "The logit-transformed propResistant intercept is -3.69"
# "The logit-transformed propResistant intercept is -3.69"

resistant_prior <- c(prior(normal(-3.7, 1), class = "Intercept"),
                     prior(normal(0, 1), class = "b"),
                     prior(normal(0, 2), class = "sds"),
                     prior(exponential(2), class = "sd"),
                     prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
                     prior(exponential(2), class = "sd", dpar = "zi"),
                     prior(normal(0, 2), class = "sds", dpar = "zi"))

resistant_prior$resp <- "propResistant"

Cyanobacteria

#-Cyanobacteria-################################################################

cyano_mod <- bf(cyanobacteria ~ 1  + rugosity_GC + mean_SR_rugosity +
                  susc_GC + mean_SR_susc_cover +
                  resistant_GC + mean_SR_resistant_cover +
                  (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                phi ~ 1  + rugosity_GC + mean_SR_rugosity +
                  susc_GC + mean_SR_susc_cover +
                  resistant_GC + mean_SR_resistant_cover +
                  (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                zi ~ 1  + rugosity_GC + mean_SR_rugosity +
                  susc_GC + mean_SR_susc_cover +
                  resistant_GC + mean_SR_resistant_cover +
                  (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

cyano_mean <- ncrmp %>%
  dplyr::group_by(SUB_REGION_NAME) %>%
  dplyr::summarize(mean = mean(cyanobacteria)) %>%
  mutate(logit_mean = logit(mean))

print(paste("The logit-transformed cyanobacteria intercept is", round(cyano_mean[1,3],2)))
## [1] "The logit-transformed cyanobacteria intercept is -1.47"
#""The logit-transformed cyanobacteria intercept is -1.47"

get_prior(cyano_mod, data = ncrmp) 
##                 prior     class                    coef           group resp
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##  student_t(3, 0, 2.5) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##  student_t(3, 0, 2.5) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##        logistic(0, 1) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##  dpar nlpar lb ub       source
##                        default
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                        default
##              0         default
##              0    (vectorized)
##              0    (vectorized)
##              0    (vectorized)
##              0    (vectorized)
##              0         default
##              0    (vectorized)
##   phi                  default
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi                  default
##   phi        0         default
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0         default
##   phi        0    (vectorized)
##    zi                  default
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi                  default
##    zi        0         default
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0         default
##    zi        0    (vectorized)
cyano_prior <- c(prior(normal(-1.5, 1), class = "Intercept"),
                      prior(normal(0, 1), class = "b"),
                      prior(exponential(2), class = "sd"),
                      prior(normal(0, 2), class = "sds"),
                      prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
                      prior(exponential(2), class = "sd", dpar = "zi"),
                      prior(normal(0, 2), class = "sds", dpar = "zi"))

cyano_prior$resp <- "cyanobacteria"   

CCA

#-CCA-##########################################################################

cca_mod <- bf(CCA ~ 1  + rugosity_GC + mean_SR_rugosity +
                susc_GC + mean_SR_susc_cover +
                resistant_GC + mean_SR_resistant_cover +
                (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
              zi ~ 1  + rugosity_GC + mean_SR_rugosity +
                susc_GC + mean_SR_susc_cover +
                resistant_GC + mean_SR_resistant_cover +
                (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
              phi ~ 1  + rugosity_GC + mean_SR_rugosity +
                susc_GC + mean_SR_susc_cover +
                resistant_GC + mean_SR_resistant_cover +
                (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
              family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

CCA_mean <- ncrmp %>%
  dplyr::group_by(SUB_REGION_NAME) %>%
  dplyr::summarize(mean = mean(CCA)) %>%
  mutate(logit_mean = logit(mean))

print(paste("The logit-transformed CCA intercept is", round(CCA_mean[1,3], 2)))
## [1] "The logit-transformed CCA intercept is -2.99"
# "The logit-transformed CCA intercept is -2.99"

get_prior(cca_mod, data = ncrmp)
##                 prior     class                    coef           group resp
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##  student_t(3, 0, 2.5) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##  student_t(3, 0, 2.5) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##        logistic(0, 1) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##  dpar nlpar lb ub       source
##                        default
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                        default
##              0         default
##              0    (vectorized)
##              0    (vectorized)
##              0    (vectorized)
##              0    (vectorized)
##              0         default
##              0    (vectorized)
##   phi                  default
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi                  default
##   phi        0         default
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0         default
##   phi        0    (vectorized)
##    zi                  default
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi                  default
##    zi        0         default
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0         default
##    zi        0    (vectorized)
cca_brm_priors <- c(prior(normal(-3, 1), class = "Intercept"),
                    prior(normal(0, 1), class = "b"),
                    prior(exponential(2), class = "sd"),
                    prior(normal(0, 2), class = "sds"),
                    prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
                    prior(exponential(2), class = "sd", dpar = "zi"),
                    prior(normal(0, 2), class = "sds", dpar = "zi"))

cca_prior <- cca_brm_priors
cca_prior$resp <- "CCA"

Macroalgae

#-Macroalgae-###################################################################

macro_mod <- bf(macroalgae ~ 1  + rugosity_GC + mean_SR_rugosity +
                  susc_GC + mean_SR_susc_cover +
                  resistant_GC + mean_SR_resistant_cover +
                  (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                zi ~ 1  + rugosity_GC + mean_SR_rugosity +
                  susc_GC + mean_SR_susc_cover +
                  resistant_GC + mean_SR_resistant_cover +
                  (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                phi ~ 1  + rugosity_GC + mean_SR_rugosity +
                  susc_GC + mean_SR_susc_cover +
                  resistant_GC + mean_SR_resistant_cover +
                  (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))

macroalgae_mean <- ncrmp %>%
  dplyr::group_by(SUB_REGION_NAME) %>%
  dplyr::summarize(mean = mean(macroalgae)) %>%
  mutate(logit_mean = logit(mean))

print(paste("The logit-transformed macroalgae intercept is", round(macroalgae_mean[1,3], 2)))
## [1] "The logit-transformed macroalgae intercept is -0.95"
# "The logit-transformed macroalgae intercept is -0.95"

get_prior(macro_mod, data = ncrmp)
##                 prior     class                    coef           group resp
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##  student_t(3, 0, 2.5) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##  student_t(3, 0, 2.5) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##        logistic(0, 1) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##  dpar nlpar lb ub       source
##                        default
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                        default
##              0         default
##              0    (vectorized)
##              0    (vectorized)
##              0    (vectorized)
##              0    (vectorized)
##              0         default
##              0    (vectorized)
##   phi                  default
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi                  default
##   phi        0         default
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0         default
##   phi        0    (vectorized)
##    zi                  default
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi                  default
##    zi        0         default
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0         default
##    zi        0    (vectorized)
macro_prior <- c(prior(normal(-1, 1), class = "Intercept"),
                      prior(normal(0, 1), class = "b"),
                      prior(exponential(2), class = "sd"),
                      prior(normal(0, 2), class = "sds"),
                      prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
                      prior(exponential(2), class = "sd", dpar = "zi"),
                      prior(normal(0, 2), class = "sds", dpar = "zi"))

macro_prior$resp <- "macroalgae"     

Fire coral

#-Fire coral-###################################################################

fire_mod <- bf(fireCoral ~ 1  + rugosity_GC + mean_SR_rugosity +
                 susc_GC + mean_SR_susc_cover +
                 resistant_GC + mean_SR_resistant_cover +
                 (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
               zi ~ 1  + rugosity_GC + mean_SR_rugosity +
                 susc_GC + mean_SR_susc_cover +
                 resistant_GC + mean_SR_resistant_cover +
                 (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
               phi ~ 1  + rugosity_GC + mean_SR_rugosity +
                 susc_GC + mean_SR_susc_cover +
                 resistant_GC + mean_SR_resistant_cover +
                 (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
               family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))


fireCoral_mean <- ncrmp %>%
  dplyr::group_by(SUB_REGION_NAME) %>%
  dplyr::summarize(mean = mean(fireCoral)) %>%
  mutate(logit_mean = logit(mean))

print(paste("The logit-transformed fireCoral intercept is", round(fireCoral_mean[1,3], 2)))
## [1] "The logit-transformed fireCoral intercept is -5.16"
# "The logit-transformed fireCoral intercept is -5.16


get_prior(fire_mod, data = ncrmp)
##                 prior     class                    coef           group resp
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##  student_t(3, 0, 2.5) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##  student_t(3, 0, 2.5) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                (flat)         b                                             
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##        logistic(0, 1) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##  dpar nlpar lb ub       source
##                        default
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                   (vectorized)
##                        default
##              0         default
##              0    (vectorized)
##              0    (vectorized)
##              0    (vectorized)
##              0    (vectorized)
##              0         default
##              0    (vectorized)
##   phi                  default
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi             (vectorized)
##   phi                  default
##   phi        0         default
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0    (vectorized)
##   phi        0         default
##   phi        0    (vectorized)
##    zi                  default
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi             (vectorized)
##    zi                  default
##    zi        0         default
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0    (vectorized)
##    zi        0         default
##    zi        0    (vectorized)
fire_prior <- c(prior(normal(-5.2, 1), class = "Intercept"),
                     prior(normal(0, 1), class = "b"),
                     prior(normal(0, 2), class = "sds"),
                     prior(exponential(2), class = "sd"),
                     prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
                     prior(exponential(2), class = "sd", dpar = "zi"),
                     prior(normal(0, 2), class = "sds", dpar = "zi"))

fire_prior$resp <- "fireCoral"

Rugosity

#-Rugosity-#####################################################################

rug_mod <- bf(logNewRug ~ 1 + susc_GC + mean_SR_susc_cover + 
                        resistant_GC + mean_SR_resistant_cover +
                        rugosityMethod + (1|SUB_REGION_NAME) + (1|YEAR),
                      sigma ~ 1 + susc_GC + mean_SR_susc_cover + 
                        resistant_GC + mean_SR_resistant_cover + 
                        rugosityMethod + (1|SUB_REGION_NAME) + (1|YEAR),family = student())


get_prior(rug_mod, data = ncrmp) 
##                  prior     class                    coef           group resp
##                 (flat)         b                                             
##                 (flat)         b mean_SR_resistant_cover                     
##                 (flat)         b      mean_SR_susc_cover                     
##                 (flat)         b            resistant_GC                     
##                 (flat)         b   rugosityMethodWTD_RUG                     
##                 (flat)         b                 susc_GC                     
##  student_t(3, -1, 2.5) Intercept                                             
##          gamma(2, 0.1)        nu                                             
##   student_t(3, 0, 2.5)        sd                                             
##   student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##   student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##   student_t(3, 0, 2.5)        sd                                    YEAR     
##   student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##                 (flat)         b                                             
##                 (flat)         b mean_SR_resistant_cover                     
##                 (flat)         b      mean_SR_susc_cover                     
##                 (flat)         b            resistant_GC                     
##                 (flat)         b   rugosityMethodWTD_RUG                     
##                 (flat)         b                 susc_GC                     
##   student_t(3, 0, 2.5) Intercept                                             
##   student_t(3, 0, 2.5)        sd                                             
##   student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##   student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##   student_t(3, 0, 2.5)        sd                                    YEAR     
##   student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##   dpar nlpar lb ub       source
##                         default
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                         default
##               1         default
##               0         default
##               0    (vectorized)
##               0    (vectorized)
##               0    (vectorized)
##               0    (vectorized)
##  sigma                  default
##  sigma             (vectorized)
##  sigma             (vectorized)
##  sigma             (vectorized)
##  sigma             (vectorized)
##  sigma             (vectorized)
##  sigma                  default
##  sigma        0         default
##  sigma        0    (vectorized)
##  sigma        0    (vectorized)
##  sigma        0    (vectorized)
##  sigma        0    (vectorized)
logNewRug_mean <- ncrmp %>%
  dplyr::group_by(SUB_REGION_NAME) %>%
  dplyr::summarize(mean = mean(logNewRug))

print(paste("The logNewRug intercept is", round(logNewRug_mean[1,2], 2)))
## [1] "The logNewRug intercept is -0.95"
# "The logNewRug intercept is -0.95

rug_prior <- c(prior(normal(-1, 1), class = "Intercept"), 
               prior(normal(0, 1), class = "b"), 
               prior(exponential(1), class = "sd"), 
               prior(normal(-1, 1), class = "Intercept", dpar = "sigma"), 
               prior(exponential(1), class = "sd", dpar = "sigma")) 

rug_prior$resp <- "logNewRug"

Coral-associated fish

#-Coral associated fish-########################################################

coralAssoc_mod<- bf(coralAssociated ~  1 + 
                      rugosity_GC+ mean_SR_rugosity + rugosityMethod +
                      susc_GC + mean_SR_susc_cover + 
                      resistant_GC + mean_SR_resistant_cover + 
                      mean_SR_cyano + cyano_GC + 
                      mean_SR_cca + cca_GC + 
                      mean_SR_macro + macro_GC + 
                      mean_SR_fire + fire_GC + 
                      (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                    shape ~ 1 + 
                      rugosity_GC+ mean_SR_rugosity + rugosityMethod +
                      susc_GC + mean_SR_susc_cover + 
                      resistant_GC + mean_SR_resistant_cover + 
                      mean_SR_cyano + cyano_GC + 
                      mean_SR_cca + cca_GC + 
                      mean_SR_macro + macro_GC + 
                      mean_SR_fire + fire_GC + 
                      (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                    hu ~ 1 + 
                      rugosity_GC+ mean_SR_rugosity + rugosityMethod +
                      susc_GC + mean_SR_susc_cover + 
                      resistant_GC + mean_SR_resistant_cover + 
                      mean_SR_cyano + cyano_GC + 
                      mean_SR_cca + cca_GC + 
                      mean_SR_macro + macro_GC + 
                      mean_SR_fire + fire_GC + 
                      (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                    family = hurdle_gamma(link = "log", link_shape = "log", link_hu = "logit"))


get_prior(coralAssoc_mod, data = ncrmp)
##                   prior     class                    coef           group resp
##                  (flat)         b                                             
##                  (flat)         b                  cca_GC                     
##                  (flat)         b                cyano_GC                     
##                  (flat)         b                 fire_GC                     
##                  (flat)         b                macro_GC                     
##                  (flat)         b             mean_SR_cca                     
##                  (flat)         b           mean_SR_cyano                     
##                  (flat)         b            mean_SR_fire                     
##                  (flat)         b           mean_SR_macro                     
##                  (flat)         b mean_SR_resistant_cover                     
##                  (flat)         b        mean_SR_rugosity                     
##                  (flat)         b      mean_SR_susc_cover                     
##                  (flat)         b            resistant_GC                     
##                  (flat)         b             rugosity_GC                     
##                  (flat)         b   rugosityMethodWTD_RUG                     
##                  (flat)         b            sMAX_DEPTH_1                     
##                  (flat)         b                 susc_GC                     
##  student_t(3, 0.4, 2.5) Intercept                                             
##    student_t(3, 0, 2.5)        sd                                             
##    student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##    student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##    student_t(3, 0, 2.5)        sd                                    YEAR     
##    student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##    student_t(3, 0, 2.5)       sds                                             
##    student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                  (flat)         b                                             
##                  (flat)         b                  cca_GC                     
##                  (flat)         b                cyano_GC                     
##                  (flat)         b                 fire_GC                     
##                  (flat)         b                macro_GC                     
##                  (flat)         b             mean_SR_cca                     
##                  (flat)         b           mean_SR_cyano                     
##                  (flat)         b            mean_SR_fire                     
##                  (flat)         b           mean_SR_macro                     
##                  (flat)         b mean_SR_resistant_cover                     
##                  (flat)         b        mean_SR_rugosity                     
##                  (flat)         b      mean_SR_susc_cover                     
##                  (flat)         b            resistant_GC                     
##                  (flat)         b             rugosity_GC                     
##                  (flat)         b   rugosityMethodWTD_RUG                     
##                  (flat)         b            sMAX_DEPTH_1                     
##                  (flat)         b                 susc_GC                     
##          logistic(0, 1) Intercept                                             
##    student_t(3, 0, 2.5)        sd                                             
##    student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##    student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##    student_t(3, 0, 2.5)        sd                                    YEAR     
##    student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##    student_t(3, 0, 2.5)       sds                                             
##    student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                  (flat)         b                                             
##                  (flat)         b                  cca_GC                     
##                  (flat)         b                cyano_GC                     
##                  (flat)         b                 fire_GC                     
##                  (flat)         b                macro_GC                     
##                  (flat)         b             mean_SR_cca                     
##                  (flat)         b           mean_SR_cyano                     
##                  (flat)         b            mean_SR_fire                     
##                  (flat)         b           mean_SR_macro                     
##                  (flat)         b mean_SR_resistant_cover                     
##                  (flat)         b        mean_SR_rugosity                     
##                  (flat)         b      mean_SR_susc_cover                     
##                  (flat)         b            resistant_GC                     
##                  (flat)         b             rugosity_GC                     
##                  (flat)         b   rugosityMethodWTD_RUG                     
##                  (flat)         b            sMAX_DEPTH_1                     
##                  (flat)         b                 susc_GC                     
##    student_t(3, 0, 2.5) Intercept                                             
##    student_t(3, 0, 2.5)        sd                                             
##    student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##    student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##    student_t(3, 0, 2.5)        sd                                    YEAR     
##    student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##    student_t(3, 0, 2.5)       sds                                             
##    student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##   dpar nlpar lb ub       source
##                         default
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                         default
##               0         default
##               0    (vectorized)
##               0    (vectorized)
##               0    (vectorized)
##               0    (vectorized)
##               0         default
##               0    (vectorized)
##     hu                  default
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu             (vectorized)
##     hu                  default
##     hu        0         default
##     hu        0    (vectorized)
##     hu        0    (vectorized)
##     hu        0    (vectorized)
##     hu        0    (vectorized)
##     hu        0         default
##     hu        0    (vectorized)
##  shape                  default
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape                  default
##  shape        0         default
##  shape        0    (vectorized)
##  shape        0    (vectorized)
##  shape        0    (vectorized)
##  shape        0    (vectorized)
##  shape        0         default
##  shape        0    (vectorized)
coralAssociated_mean <- ncrmp %>%
  dplyr::group_by(SUB_REGION_NAME) %>%
  dplyr::summarize(mean = mean(coralAssociated)) %>%
  mutate(logmean = log(mean))

print(paste("The log-transformed coralAssociated intercept is", round(coralAssociated_mean[1,3], 2)))
## [1] "The log-transformed coralAssociated intercept is 2.31"
# "The log-transformed coralAssociated intercept is 2.31"


coralAssoc_prior <- c(prior(normal(2.3, 1), class = "Intercept"),
                      prior(normal(0, 1), class = "b"),
                      prior(normal(0, 2), class = "sds"),
                      prior(exponential(1), class = "sd"),
                      prior(normal(2.3, 1), class = "Intercept", dpar = "shape"),
                      prior(exponential(1), class = "sd", dpar = "shape"), 
                      prior(normal(0, 1), class = "b", dpar = "shape"))
coralAssoc_prior$resp <- "coralAssociated"

Herbivorous fish

#-Herbivorous fishes-###########################################################

herb_mod <- bf(herb ~ 1 + mean_SR_rugosity + rugosity_GC + 
                 susc_GC + mean_SR_susc_cover + 
                 resistant_GC + mean_SR_resistant_cover + 
                 mean_SR_cyano + cyano_GC + 
                 mean_SR_cca + cca_GC + 
                 mean_SR_macro + macro_GC + 
                 mean_SR_fire + fire_GC + 
                 (1 | SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
               shape ~ 1 + mean_SR_rugosity + rugosity_GC + 
                 susc_GC + mean_SR_susc_cover + 
                 resistant_GC + mean_SR_resistant_cover + 
                 mean_SR_cyano + cyano_GC + 
                 mean_SR_cca + cca_GC + 
                 mean_SR_macro + macro_GC + 
                 mean_SR_fire + fire_GC + 
                 (1 | SUB_REGION_NAME) + (1|YEAR)+ s(MAX_DEPTH),   
               family = hurdle_gamma(link = "log"))


get_prior(herb_mod, data = ncrmp)
##                    prior     class                    coef           group resp
##                   (flat)         b                                             
##                   (flat)         b                  cca_GC                     
##                   (flat)         b                cyano_GC                     
##                   (flat)         b                 fire_GC                     
##                   (flat)         b                macro_GC                     
##                   (flat)         b             mean_SR_cca                     
##                   (flat)         b           mean_SR_cyano                     
##                   (flat)         b            mean_SR_fire                     
##                   (flat)         b           mean_SR_macro                     
##                   (flat)         b mean_SR_resistant_cover                     
##                   (flat)         b        mean_SR_rugosity                     
##                   (flat)         b      mean_SR_susc_cover                     
##                   (flat)         b            resistant_GC                     
##                   (flat)         b             rugosity_GC                     
##                   (flat)         b            sMAX_DEPTH_1                     
##                   (flat)         b                 susc_GC                     
##               beta(1, 1)        hu                                             
##  student_t(3, -0.6, 2.5) Intercept                                             
##     student_t(3, 0, 2.5)        sd                                             
##     student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##     student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##     student_t(3, 0, 2.5)        sd                                    YEAR     
##     student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##     student_t(3, 0, 2.5)       sds                                             
##     student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                   (flat)         b                                             
##                   (flat)         b                  cca_GC                     
##                   (flat)         b                cyano_GC                     
##                   (flat)         b                 fire_GC                     
##                   (flat)         b                macro_GC                     
##                   (flat)         b             mean_SR_cca                     
##                   (flat)         b           mean_SR_cyano                     
##                   (flat)         b            mean_SR_fire                     
##                   (flat)         b           mean_SR_macro                     
##                   (flat)         b mean_SR_resistant_cover                     
##                   (flat)         b        mean_SR_rugosity                     
##                   (flat)         b      mean_SR_susc_cover                     
##                   (flat)         b            resistant_GC                     
##                   (flat)         b             rugosity_GC                     
##                   (flat)         b            sMAX_DEPTH_1                     
##                   (flat)         b                 susc_GC                     
##     student_t(3, 0, 2.5) Intercept                                             
##     student_t(3, 0, 2.5)        sd                                             
##     student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##     student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##     student_t(3, 0, 2.5)        sd                                    YEAR     
##     student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##     student_t(3, 0, 2.5)       sds                                             
##     student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##   dpar nlpar lb ub       source
##                         default
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##               0  1      default
##                         default
##               0         default
##               0    (vectorized)
##               0    (vectorized)
##               0    (vectorized)
##               0    (vectorized)
##               0         default
##               0    (vectorized)
##  shape                  default
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape                  default
##  shape        0         default
##  shape        0    (vectorized)
##  shape        0    (vectorized)
##  shape        0    (vectorized)
##  shape        0    (vectorized)
##  shape        0         default
##  shape        0    (vectorized)
herb_mean <- ncrmp %>%
  dplyr::group_by(SUB_REGION_NAME) %>%
  dplyr::summarize(mean = mean(herb)) %>%
  mutate(logmean = log(mean))

print(paste("The log-transformed herb intercept is", round(herb_mean[1,3], 2)))
## [1] "The log-transformed herb intercept is -0.39"
# "The log-transformed herb intercept is -0.39"

herb_prior <- c(prior(normal(-0.4, 1), class = "Intercept"),
                prior(normal(0, 1), class = "b"),
                prior(normal(0, 2), class = "sds"),
                prior(exponential(1), class = "sd"),
                prior(beta(1, 1), class = "hu"),
                prior(normal(-0.4, 1), class = "Intercept", dpar = "shape"),
                prior(exponential(1), class = "sd", dpar = "shape"), 
                prior(normal(0, 1), class = "b", dpar = "shape"))


herb_prior$resp <- "herb"

Other fish

#-Other fish-###################################################################

otherFish_mod <- bf(otherFish ~ 1 + 
                            rugosity_GC+ mean_SR_rugosity + rugosityMethod +
                            susc_GC + mean_SR_susc_cover + 
                            resistant_GC + mean_SR_resistant_cover + 
                            mean_SR_cyano + cyano_GC + 
                            mean_SR_cca + cca_GC + 
                            mean_SR_macro + macro_GC + 
                            mean_SR_fire + fire_GC + 
                            (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                          shape ~ 1 + 
                            rugosity_GC+ mean_SR_rugosity + rugosityMethod +
                            susc_GC + mean_SR_susc_cover + 
                            resistant_GC + mean_SR_resistant_cover + 
                            mean_SR_cyano + cyano_GC + 
                            mean_SR_cca + cca_GC + 
                            mean_SR_macro + macro_GC + 
                            mean_SR_fire + fire_GC + 
                            (1|SUB_REGION_NAME) + (1|YEAR) + s(MAX_DEPTH),
                          family = Gamma(link = "log"))

get_prior(otherFish_mod, data = ncrmp)
##                 prior     class                    coef           group resp
##                (flat)         b                                             
##                (flat)         b                  cca_GC                     
##                (flat)         b                cyano_GC                     
##                (flat)         b                 fire_GC                     
##                (flat)         b                macro_GC                     
##                (flat)         b             mean_SR_cca                     
##                (flat)         b           mean_SR_cyano                     
##                (flat)         b            mean_SR_fire                     
##                (flat)         b           mean_SR_macro                     
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b   rugosityMethodWTD_RUG                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##  student_t(3, 1, 2.5) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##                (flat)         b                                             
##                (flat)         b                  cca_GC                     
##                (flat)         b                cyano_GC                     
##                (flat)         b                 fire_GC                     
##                (flat)         b                macro_GC                     
##                (flat)         b             mean_SR_cca                     
##                (flat)         b           mean_SR_cyano                     
##                (flat)         b            mean_SR_fire                     
##                (flat)         b           mean_SR_macro                     
##                (flat)         b mean_SR_resistant_cover                     
##                (flat)         b        mean_SR_rugosity                     
##                (flat)         b      mean_SR_susc_cover                     
##                (flat)         b            resistant_GC                     
##                (flat)         b             rugosity_GC                     
##                (flat)         b   rugosityMethodWTD_RUG                     
##                (flat)         b            sMAX_DEPTH_1                     
##                (flat)         b                 susc_GC                     
##  student_t(3, 0, 2.5) Intercept                                             
##  student_t(3, 0, 2.5)        sd                                             
##  student_t(3, 0, 2.5)        sd                         SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd               Intercept SUB_REGION_NAME     
##  student_t(3, 0, 2.5)        sd                                    YEAR     
##  student_t(3, 0, 2.5)        sd               Intercept            YEAR     
##  student_t(3, 0, 2.5)       sds                                             
##  student_t(3, 0, 2.5)       sds            s(MAX_DEPTH)                     
##   dpar nlpar lb ub       source
##                         default
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                    (vectorized)
##                         default
##               0         default
##               0    (vectorized)
##               0    (vectorized)
##               0    (vectorized)
##               0    (vectorized)
##               0         default
##               0    (vectorized)
##  shape                  default
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape             (vectorized)
##  shape                  default
##  shape        0         default
##  shape        0    (vectorized)
##  shape        0    (vectorized)
##  shape        0    (vectorized)
##  shape        0    (vectorized)
##  shape        0         default
##  shape        0    (vectorized)
otherFish_mean <- ncrmp %>%
  dplyr::group_by(SUB_REGION_NAME) %>%
  dplyr::summarize(mean = mean(otherFish)) %>%
  mutate(logmean = log(mean))

print(paste("The log-transformed otherFish intercept is", round(otherFish_mean[1,3], 2)))
## [1] "The log-transformed otherFish intercept is 3.49"
# "The log-transformed otherFish intercept is 3.49

otherFish_prior <- c(prior(normal(3.5, 1), class = "Intercept"),
                     prior(normal(0, 1), class = "b"),
                     prior(normal(0, 2), class = "sds"),
                     prior(exponential(1), class = "sd"),
                     prior(normal(3.1, 1), class = "Intercept", dpar ="shape"), 
                     prior(normal(0, 1), class = "b", dpar = "shape"),
                     prior(exponential(1), class = "sd", dpar = "shape"))

otherFish_prior$resp <- "otherFish"

Build and run full SEM

################################################################################
#-Build full SEM-###############################################################
################################################################################

full_priors <- c( susc_prior,
                  resistant_prior,
                  cyano_prior, 
                 cca_prior,
                 macro_prior, 
                 fire_prior, 
                 rug_prior, 
                 coralAssoc_prior,
                herb_prior,
                otherFish_prior)
full_priors
##            prior     class coef group            resp  dpar nlpar   lb   ub
##  normal(-2.1, 1) Intercept                   propSusc             <NA> <NA>
##     normal(0, 1)         b                   propSusc             <NA> <NA>
##     normal(0, 2)       sds                   propSusc             <NA> <NA>
##   exponential(2)        sd                   propSusc             <NA> <NA>
##   logistic(0, 1) Intercept                   propSusc    zi       <NA> <NA>
##   exponential(2)        sd                   propSusc    zi       <NA> <NA>
##     normal(0, 2)       sds                   propSusc    zi       <NA> <NA>
##  normal(-3.7, 1) Intercept              propResistant             <NA> <NA>
##     normal(0, 1)         b              propResistant             <NA> <NA>
##     normal(0, 2)       sds              propResistant             <NA> <NA>
##   exponential(2)        sd              propResistant             <NA> <NA>
##   logistic(0, 1) Intercept              propResistant    zi       <NA> <NA>
##   exponential(2)        sd              propResistant    zi       <NA> <NA>
##     normal(0, 2)       sds              propResistant    zi       <NA> <NA>
##  normal(-1.5, 1) Intercept              cyanobacteria             <NA> <NA>
##     normal(0, 1)         b              cyanobacteria             <NA> <NA>
##   exponential(2)        sd              cyanobacteria             <NA> <NA>
##     normal(0, 2)       sds              cyanobacteria             <NA> <NA>
##   logistic(0, 1) Intercept              cyanobacteria    zi       <NA> <NA>
##   exponential(2)        sd              cyanobacteria    zi       <NA> <NA>
##     normal(0, 2)       sds              cyanobacteria    zi       <NA> <NA>
##    normal(-3, 1) Intercept                        CCA             <NA> <NA>
##     normal(0, 1)         b                        CCA             <NA> <NA>
##   exponential(2)        sd                        CCA             <NA> <NA>
##     normal(0, 2)       sds                        CCA             <NA> <NA>
##   logistic(0, 1) Intercept                        CCA    zi       <NA> <NA>
##   exponential(2)        sd                        CCA    zi       <NA> <NA>
##     normal(0, 2)       sds                        CCA    zi       <NA> <NA>
##    normal(-1, 1) Intercept                 macroalgae             <NA> <NA>
##     normal(0, 1)         b                 macroalgae             <NA> <NA>
##   exponential(2)        sd                 macroalgae             <NA> <NA>
##     normal(0, 2)       sds                 macroalgae             <NA> <NA>
##   logistic(0, 1) Intercept                 macroalgae    zi       <NA> <NA>
##   exponential(2)        sd                 macroalgae    zi       <NA> <NA>
##     normal(0, 2)       sds                 macroalgae    zi       <NA> <NA>
##  normal(-5.2, 1) Intercept                  fireCoral             <NA> <NA>
##     normal(0, 1)         b                  fireCoral             <NA> <NA>
##     normal(0, 2)       sds                  fireCoral             <NA> <NA>
##   exponential(2)        sd                  fireCoral             <NA> <NA>
##   logistic(0, 1) Intercept                  fireCoral    zi       <NA> <NA>
##   exponential(2)        sd                  fireCoral    zi       <NA> <NA>
##     normal(0, 2)       sds                  fireCoral    zi       <NA> <NA>
##    normal(-1, 1) Intercept                  logNewRug             <NA> <NA>
##     normal(0, 1)         b                  logNewRug             <NA> <NA>
##   exponential(1)        sd                  logNewRug             <NA> <NA>
##    normal(-1, 1) Intercept                  logNewRug sigma       <NA> <NA>
##   exponential(1)        sd                  logNewRug sigma       <NA> <NA>
##   normal(2.3, 1) Intercept            coralAssociated             <NA> <NA>
##     normal(0, 1)         b            coralAssociated             <NA> <NA>
##     normal(0, 2)       sds            coralAssociated             <NA> <NA>
##   exponential(1)        sd            coralAssociated             <NA> <NA>
##   normal(2.3, 1) Intercept            coralAssociated shape       <NA> <NA>
##   exponential(1)        sd            coralAssociated shape       <NA> <NA>
##     normal(0, 1)         b            coralAssociated shape       <NA> <NA>
##  normal(-0.4, 1) Intercept                       herb             <NA> <NA>
##     normal(0, 1)         b                       herb             <NA> <NA>
##     normal(0, 2)       sds                       herb             <NA> <NA>
##   exponential(1)        sd                       herb             <NA> <NA>
##       beta(1, 1)        hu                       herb             <NA> <NA>
##  normal(-0.4, 1) Intercept                       herb shape       <NA> <NA>
##   exponential(1)        sd                       herb shape       <NA> <NA>
##     normal(0, 1)         b                       herb shape       <NA> <NA>
##   normal(3.5, 1) Intercept                  otherFish             <NA> <NA>
##     normal(0, 1)         b                  otherFish             <NA> <NA>
##     normal(0, 2)       sds                  otherFish             <NA> <NA>
##   exponential(1)        sd                  otherFish             <NA> <NA>
##   normal(3.1, 1) Intercept                  otherFish shape       <NA> <NA>
##     normal(0, 1)         b                  otherFish shape       <NA> <NA>
##   exponential(1)        sd                  otherFish shape       <NA> <NA>
##  source
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
##    user
### UNCOMMENT BELOW TO RUN MODEL
# full_SEM <- brm(cyano_mod + 
#                 rug_mod + 
#                 cca_mod + 
#                 fire_mod + 
#                 macro_mod + 
#                 otherFish_mod + 
#                 coralAssoc_mod + 
#                 herb_mod +
#                 resistantCoral_mod + 
#                 suscCoral_mod +
#                 set_rescor(rescor = FALSE),
#                         prior = full_priors,
#                         data = ncrmp,
#                         cores = n_cores,
#                         chains = n_chains,
#                         iter = n_iter,
#                         warmup = n_warmup,
#                         init = "0",
#                         control = list(adapt_delta = 0.99))

## save to file
full_SEM <- readRDS("data/outputs/full_SEM.rds")
#saveRDS(full_SEM, "data/outputs/full_SEM.rds")

Posterior predictive checks

SCTLD-susceptible coral

pp_check(full_SEM, ndraws = 100,resp = 'propSusc', type = "dens_overlay")+ 
  labs(title = "Sctld-susceptible coral")

pp_check(full_SEM, ndraws = 100,resp = 'propSusc', type = "hist")+ 
  labs(title = "Sctld-susceptible coral")

pp_check(full_SEM, ndraws = 100,resp = 'propSusc', type = "boxplot")+ 
  labs(title = "Sctld-susceptible coral")

y <- t(posterior_predict(full_SEM, resp = 'propSusc')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$propSusc*100)
length(y_int) 
## [1] 1637
ncol(yrep_int)
## [1] 1637
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

SCTLD-resistant coral

pp_check(full_SEM, ndraws = 100,resp = 'propResistant', type = "dens_overlay")+ 
  labs(title = "Sctld-resistant coral")

pp_check(full_SEM, ndraws = 100,resp = 'propResistant', type = "hist")+ 
  labs(title = "Sctld-resistant coral")

pp_check(full_SEM, ndraws = 100,resp = 'propResistant', type = "boxplot")+ 
  labs(title = "Sctld-resistant coral")

y <- t(posterior_predict(full_SEM, resp = 'propResistant')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$propResistant*100)
length(y_int) 
## [1] 1637
ncol(yrep_int)
## [1] 1637
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

cyanobacteria

pp_check(full_SEM, ndraws = 100,resp = 'cyanobacteria', type = "dens_overlay")+ 
  labs(title = "cyanobacteria")

pp_check(full_SEM, ndraws = 100,resp = 'cyanobacteria', type = "hist")+ 
  labs(title = "cyanobacteria")

pp_check(full_SEM, ndraws = 100,resp = 'cyanobacteria', type = "boxplot")+ 
  labs(title = "cyanobacteria")

y <- t(posterior_predict(full_SEM, resp = 'cyanobacteria')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$cyanobacteria*100)
length(y_int) 
## [1] 1637
ncol(yrep_int)
## [1] 1637
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

CCA

pp_check(full_SEM, ndraws = 100,resp = 'CCA', type = "dens_overlay")+ 
  labs(title = "cca")

pp_check(full_SEM, ndraws = 100,resp = 'CCA', type = "hist")+ 
  labs(title = "cca")

pp_check(full_SEM, ndraws = 100,resp = 'CCA', type = "boxplot")+ 
  labs(title = "cca")

y <- t(posterior_predict(full_SEM, resp = 'CCA')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$CCA*100)
length(y_int) 
## [1] 1637
ncol(yrep_int)
## [1] 1637
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

macroalgae

pp_check(full_SEM, ndraws = 100,resp = 'macroalgae', type = "dens_overlay")+ 
  labs(title = "macroalgae")

pp_check(full_SEM, ndraws = 100,resp = 'macroalgae', type = "hist")+ 
  labs(title = "macroalgae")

pp_check(full_SEM, ndraws = 100,resp = 'macroalgae', type = "boxplot")+ 
  labs(title = "macroalgae")

y <- t(posterior_predict(full_SEM, resp = 'macroalgae')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$macroalgae*100)
length(y_int) 
## [1] 1637
ncol(yrep_int)
## [1] 1637
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

fire coral

pp_check(full_SEM, ndraws = 100,resp = 'fireCoral', type = "dens_overlay")+ 
 
  labs(title = "fire coral")

pp_check(full_SEM, ndraws = 100,resp = 'fireCoral', type = "hist")+ 
 
  labs(title = "fire coral")

pp_check(full_SEM, ndraws = 100,resp = 'fireCoral', type = "boxplot")+ 
 
  labs(title = "fire coral")

y <- t(posterior_predict(full_SEM, resp = 'fireCoral')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$fireCoral*100)
length(y_int) 
## [1] 1637
ncol(yrep_int)
## [1] 1637
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

rugosity

pp_check(full_SEM, ndraws = 100,resp = 'logNewRug', type = "dens_overlay")+ 
  labs(title = "rugosity")

pp_check(full_SEM, ndraws = 100,resp = 'logNewRug', type = "hist")+ 
  labs(title = "rugosity")

pp_check(full_SEM, ndraws = 100,resp = 'logNewRug', type = "boxplot")+ 
  labs(title = "rugosity")

pp_check(full_SEM, ndraws = 100,resp = 'logNewRug', type = "loo_pit_overlay")+ 
  labs(title = "rugosity")
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

coral-associated fish

pp_check(full_SEM, ndraws = 100,resp = 'coralAssociated', type = "dens_overlay")+ 
  labs(title = "coral-associated fish")

pp_check(full_SEM, ndraws = 100,resp = 'coralAssociated', type = "hist")+ 
  labs(title = "coral-associated fish")

pp_check(full_SEM, ndraws = 100,resp = 'coralAssociated', type = "boxplot")+ 
  labs(title = "coral-associated fish")

y <- t(posterior_predict(full_SEM, resp = 'coralAssociated')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$fireCoral*100)
length(y_int) 
## [1] 1637
ncol(yrep_int)
## [1] 1637
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

herbivorous fish

pp_check(full_SEM, ndraws = 100,resp = 'herb', type = "dens_overlay")+ 
  labs(title = "herbivorous fish")

pp_check(full_SEM, ndraws = 100,resp = 'herb', type = "hist")+ 
  labs(title = "herbivorous fish")

pp_check(full_SEM, ndraws = 100,resp = 'herb', type = "boxplot")+ 
  labs(title = "herbivorous fish")

y <- t(posterior_predict(full_SEM, resp = 'herb')*100)
all_col_indices <- 1:ncol(y)
# Randomly select 1637 row indices from the vector of all indices
selected_col_indices <- sample(all_col_indices, size = 1637)
# Subset the DataFrame using the selected row indices
selected_cols <- y[,selected_col_indices ]
yrep_int <- sapply(data.frame(selected_cols), as.integer)
y_int <- as.integer(ncrmp$fireCoral*100)
length(y_int) 
## [1] 1637
ncol(yrep_int)
## [1] 1637
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

other fish

pp_check(full_SEM, ndraws = 100,resp = 'otherFish', type = "dens_overlay")+ 
  labs(title = "otherFish")

pp_check(full_SEM, ndraws = 100,resp = 'otherFish', type = "hist")+ 
  labs(title = "otherFish")

pp_check(full_SEM, ndraws = 100,resp = 'otherFish', type = "boxplot")+ 
  labs(title = "otherFish")

pp_check(full_SEM, ndraws = 100,resp = 'otherFish', type = "loo_pit_overlay")+ 
  labs(title = "otherFish")
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.